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Oh Abstract 

2 We estimate the Lyapunov times (characteristic times of predictabil- 

ity ity of motion) in Quillen's models for the dynamics in the solar neigh- 

i borhood. These models take into account perturbations due to the 

Galactic bar and spiral arms. For estimating the Lyapunov times, an 
approach based on the separatrix map theory is used. The Lyapunov 
times turn out to be typically of the order of 10 Galactic years. We 
show that only in a narrow range of possible values of the problem 
^ parameters the Galactic chaos is adiabatic; usually it is not slow. We 

also estimate the characteristic diffusion times in the chaotic domain. 
In a number of models, the diffusion times turn out to be small enough 
to permit migration of the Sun from the inner regions of the Milky 
Way to its current location. Moreover, due to the possibility of ballis- 
tic flights inside the chaotic layer, the chaotic mixing might be even far 
^ more effective and quicker than in the case of normal diffusion. This 

^ confirms the dynamical possibility of Minchev and Famaey's migration 

concept. 
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1 Introduction 



Chaotic dynamics due to interaction of nonlinear resonances in Hamiltonian 
systems is studied in a broad range of application areas, from plasma physics 
to celestial mechanics 



(e.g. 



Chirikov 


1979; 


Lichtenberg & Lieberman 


1992 



Abdullaev 2006). The characteristic time of predictability of any motion is 



nothing but the Lyapunov time (the inverse of the maximum Lyapunov expo- 
nent) of the motion. Generally, the estimation of the Lyapunov exponents is 



one of the most important tools in the study of chaotic motion (Lichtenberg 



& Lieberman, 1992), in particular in celestial mechanics. The Lyapunov ex- 



ponents characterize the mean rate of exponential divergence of trajectories 
close to each other in phase space. Non-zero Lyapunov exponents indicate 
chaotic character of motion, while the maximum Lyapunov exponent equal 
to zero signifies regular (periodic or quasiperiodic) motion. 

The development of methods of numerical computation of the Lyapunov 
exponents has more than a 30 year history (e.g., Froeschle||1984 Lichtenberg 



&: Lieberman|1992 ). In contrast, methods of analytical estimation of the Lya- 



punov exponents started to be developed only recently (Holman & Murray 



1996 Murray & Holman 1997; Shevchenko, 2000a 20021 2007 2008b) 



In studies of the dynamics of the Milky Way, analysis of the Lyapunov 
exponents has not yet been widely used, but nevertheless there are important 
achievements. Fux (2001) used Lyapunov exponents as a tool to find bar- 
induced manifestations of chaos in the local disk stellar kinematics. Taking 
into account the perturbations from the bar solely (for some particular bar 
strengths), he constructed "Lyapunov diagrams", presenting the Lyapunov 
timescales of the orbits in the u - v velocity plane at fixed space positions, 
and identified regular and chaotic domains in velocity space as a function 
of space position with respect to the bar. The Lyapunov exponents were 
calculated on a Cartesian grid of planar velocities. The fraction of chaotic 
orbits was demonstrated to obviously increase with the bar strength. How- 
ever, the diagrams in Fux (2001 ) can hardly be used to estimate real values of 
Lyapunov times, because the saturation of the computed values of the Lya- 
punov exponents was not controlled, while the adopted computation time, 
corresponding to three Hubble times (equivalently, ~ 100 Galactic years), 
might not be enough for the saturation. In other words, the obtained val- 
ues characterize not the whole chaotic regions of the phase space, but only 
rather small vicinities of the initial data. Therefore, the computed values 
are the "local" values of the Lyapunov exponents. This is testified by the 
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fact that the variation of the Lyapunov exponents in the diagrams of Fux 



( 2001[ ) is smooth, while there must be a sharp distinction between the chaotic 
regions (with non-zero Lyapunov exponents) and regular regions (with zero 
Lyapunov exponents) in the divided phase space. 

In connection with the problem of estimation of Lyapunov timescales 



in the solar neighborhood, Quillen (2003) noted that, according to Holman 



& Murray (1996), for a fully overlapped system, the chaotic zone should 
have a Lyapunov time ~ 2ir/v (where v is the frequency of perturbation), 
corresponding to the separatrix pulsation period. In what follows we shall 
consider the model set of Quillen (2003) for the stellar dynamics in the solar 
neighborhood and show that the heuristic estimate ~ 2ir/v severely underes- 
timates the real Lyapunov time. Besides, we shall see that it is rather seldom 
that the considered dynamical systems, modeling the dynamics in the solar 
neighborhood, can be called "fully overlapped" . 

Besides obtaining the Lyapunov times, we shall estimate the diffusion 
times in the chaotic domain of the phase space, in the same model set. This 
will allow one to judge on the possibility for migration of the Sun from the 
inner regions of the Milky Way to its current location. Such a possibility, 
arising due to the overlapping of resonances in the phase space, was advocated 
and studied in detail by Minchev & Famaey (2010) and Minchev et al. (2011 ). 



2 The model of interaction of nonlinear res- 
onances 

Many problems on nonlinear resonances in mechanics and physics are de- 
scribed by the perturbed pendulum Hamiltonian 

Qp 2 

H = — J 7 cos + a cos(p — t) + b cos(ip + t) (1) 



e.g., Shevchenko 2000b). The first two terms in equation rtTj) represent the 



Hamiltonian H of the unperturbed pendulum; (p is the pendulum angle (the 
resonance phase angle) and p is the momentum. The periodic perturbations 
are given by the last two terms; r is the phase angle of perturbation: r = 
Qt + r , where Q is the perturbation frequency and r is the initial phase of 
the perturbation. The quantities J 7 , Q, a, b are parameters. 

Generally, equation ([I]) describes a triplet (triad) of resonances: there are 
three trigonometric terms, each corresponding to a particular resonance. In 
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the following Sections, the case of a resonance duad is considered (i.e., a or 
b equals zero). 

It is convenient to describe the motion in the vicinity of the separatrices 
of the Hamiltonian (II]) by means of the so-called separatrix map ( Chirikov 
1979). It is a two-dimensional, area-preserving map given by 



w i+1 = Wi- W sin n, 
32 

Ti+i = Tj + A In | 1 (mod 2ir). 



(2) 



These equations give the classical separatrix map (Chirikov, 1979), valid in 
the symmetric (a = b) case of perturbation. The variable w of the map 
denotes the relative (with respect to the unperturbed separatrix value) pen- 
dulum energy w = ^ — 1, and r retains its meaning of the phase angle 
of perturbation. The parameter A is the ratio of the absolute value of 
the perturbation frequency, to tu = \J r Q\ 1 ^ 2 , the frequency of the small- 
amplitude pendulum oscillations. The parameter W in the case of a = b has 



the form (Shevchenko, 1998) 



W = e\ (A 2 (\) + A 2 {-\)) 



AiteX 2 
sinh ^ ' 



(3) 



where e = a/ J 7 , and A 2 (A) is the value of the Melnikov- Arnold integral as 



defined in Chirikov (1979): 



A 2 (A) = 4vrA 



exp(7rA/2) 



(4) 



sinh(7rA) 

Formula (|3~] ) differs from that given in Chirikov (1979) and Lichtenberg & 



Lieberman (1992) by the term A 2 (— A), which is small for A ^> 1. However, 



its contribution is significant for small values of A (Shevchenko, 1998), i.e., in 



the case of adiabatic chaos. The kind of chaos (adiabatic or non-adiabatic) 
in model ([T]) is identified by the value of A: if A < 1/2, it is slow (adiabatic), 
otherwise it is "fast" (non-adiabatic; Shevchenko 2008a). 



One iteration of the separatrix map corresponds to one period of the 
pendulum rotation or a half-period of its libration. The applicability of the 
theory of separatrix maps for description of the motion near the separatrices 
of the perturbed nonlinear resonance in the full range of the relative frequency 
of perturbation, including its low values, was considered and shown to be 



legitimate in Shevchenko (2000b). 
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The separatrix map in the case of asymmetric (a 7^ b) perturbation is 
different from that in the symmetric (a = b) case, because the energy in- 
crements are different for the prograde and the retrograde motions of the 



model pendulum (Shevchenko, 1999). (The motion is called prograde or ret- 



rograde if the variation of (p with time is, respectively, positive or negative.) 
The algorithm, taking this difference into account, constitutes the separatrix 



algorithmic map (Shevchenko, 1999): 



if Wi < and W = W 
if Wi < and W = W 



W i+ i = Wi 



n+1 = n + A In 



32 



\w 



then W : = 
then W : = 

(mod 27r), 



w 



i+ll 



(5) 



where W + and W~ denote the values of the W parameter for the prograde 
and retrograde motions respectively. In the case of asymmetric perturbation 
these values are different. 

Equations ([5]) as well can be written in a shorter way (Shevchenko, 2000b): 



if Wi < and W = W ± then W := W T , 

w i+ i =Wi- W sin Ti, 
32 

r i+ i = Ti + A In 1 (mod 2ir). 



\w 



(6) 



i+ll 



The sign in the upper index of W alternates at each iteration if Wi < (i.e., 
at librations); the quantity W ± denotes W + or W~, while W T denotes a 
corresponding value of W~ or W + . 

The essence of the separatrix algorithmic map is in taking into account the 
alternations of the W parameter. It alternates when the direction of motion 
alternates. This takes place either when rotation changes to libration, or 
when the motion is librational. Algorithms (|5| and do not contain the 
condition w,i > 0, because the direction of motion does not change when it 
holds. 

In order to find expressions for W + and W~ , one should integrate the 
increment of energy per one iteration of the map, following the usual proce- 



dure (Chirikov, 1979), but making it separately for prograde and retrograde 



directions of motion. This gives (Shevchenko, 1999) 
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W + (X,r ] )=eX(A 2 (X) + r ] A 2 (-X)), (7) 

W~(X, rj) = eX (r)Az(X) + A 2 {-X)) , (8) 

where e = a/J r and rj = b/a. The Melnikov- Arnold integral ^(A) is given 
by equation ([5]). 

The separatrix map theory can be used for analytical estimation of the 



Shevchenko 


2000a 


2002 


2007 


2008b) 



parisons of the predictions of this theory with numerical-experimental results 



can be found in Shevchenko & Kouprianov (2002) and Shevchenko (2000a 



2002, 2007, 2008b 2009), where it was applied to various problems of celes- 
tial mechanics: rotational dynamics of planetary satellites, orbital dynamics 
of satellite systems, and orbital dynamics of asteroids. 

The maximum Lyapunov exponent is defined by the limit 



L = lim sup 



1 



t-tn 



In 



d(t )^0 



d(t) 
d(t ) 



(9) 



where c?(t ) is the distance (in phase space) between two nearby initial condi- 
tions for two trajectories at the initial moment of time to, d(t) is the distance 
between the evolved initial conditions at time t (e.g., Lichtenberg & Lieber- 



man 



1992). 



According to the general approach used by Shevchenko (2000a 
the maximum Lyapunov exponent 



2002 



mam 



2007), the maximum Lyapunov exponent, L, of the motion in the 
chaotic layer of system given by Eq. ([T| is represented as the ratio of the 
maximum Lyapunov exponent, L sx , of its separatrix map and the average 
period T of rotation (or, equivalently, the average half-period of libration) 
of the resonance phase angle <p inside the layer: L = L sx /T. In this way, 



formulas for the Lyapunov time were derived in Shevchenko (2007) for four 
generic cases of interacting resonances: the fastly chaotic resonance triplet, 
fastly chaotic resonance doublet, slowly chaotic resonance triplet, and slowly 
chaotic resonance doublet (called, respectively, "ft," "fd," "st," and "sd" 
resonance multiplet types). 

In the following Sections, we shall need formulas for the Lyapunov times 
for the "fd" and "sd" resonance types only. The formula for the Lyapunov 
time for the "fd" resonance type is 
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pert 
~2tT 



Mlibr + 1 

(2A) , 



T ax (2A,iy) + T SX (X,W) 



(10) 



= 27r/|fi| is the 



where /iubr ~ 4 (Shevchenko, 2007). The quantity T pcrt 
period of perturbation. The quantities W, L sx , and T sx are given by the 
formulas 



W{e,X) = eX (A 2 (X) + A 2 (-X)) 
L SX (X) ~ Ch 



4neX 2 



sinh - c 



A > 



T sx (A,W/)^Aln 



1 + 2A' 
32e 



Ai^r 



(ii) 

(12) 
(13) 



where Ch ~ 0.80 is Chirikov's constant (Shevchenko, 2004) and e is the base 
of natural logarithms. 



The formula for the Lyapunov time for the "sd" resonance type is 



pert 

~2^T 



In 



32 

1x 



sin 




(14) 



(Shevchenko, 2007). 



The parameter A = |fi|/u;o measures the separation of the perturbing and 
guiding resonances in the units of one quarter of the width of the guiding 
resonance, because the resonances separation in frequency space is equal to 



f2, while the guiding resonance width is equal to Aujq (Chirikov, 1979). There 



fore, A can be regarded as a kind of the resonance overlap parameter. In the 
asymptotic limit of the adiabatic (slow) case the resonances in the multiplet 
strongly overlap, while in the asymptotic limit of the "fast" case the reso- 
nances are segregated. However note that the border A ~ 1/2 (Shevchenko 



2008a) between slow and fast chaos does not coincide with the borderline 



between the cases of overlapping and non-overlapping of resonances: the lat- 
ter borderline lies much higher in A. For example, in the phase space of the 
standard map the integer resonances start to overlap (on decreasing A) at 
A w 2tt/0.97 w 6.5 (|Chirikov[ [l979|. 
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3 The resonance Hamiltonian 



Quillen (2003, 2009), based on the dynamical model of Contopoulos (1975 



1988), constructed a Hamiltonian of the motion in the solar neighborhood by 



adding the perturbations due to the bar and spiral arms to the unperturbed 
Hamiltonian of the motion. Namely, Quillen's model describes interaction 
of the bar's 2:1 outer Lindblad resonance with the spiral's 2:1 or 4:1 inner 
Lindblad resonance. The resulting Hamiltonian has the form 



H = A[f + Sj + fij 



1/2 



COS < 



+ ej 1/2 cos( 



+ vt- 7)] 



(15) 



(Quillen 2003, equation (23)). Here, j and 4> are the conjugate momentum 
and phase variables; A rs —5.7; 5, (3, e, v, and 7 are free unitless parameters. 



The ranges for numerical values of the parameters were estimated in Quillen 



(2003) from observational physical and kinematical considerations. 

The frequency v is counted in units of Q, and, accordingly, time is in 
units of f2 _1 ; here Q is the rotation rate of the epicyclic center. Therefore, 
one time unit is equal to one Galactic year at a given distance from the center 
of the Milky Way, divided by 2n. 



The first resonant term (that with the coefficient P) in equation (15) cor- 



responds to the perturbation due to the bar, while the second one (that with 
the coefficient e) to the perturbation due to the spiral arms. The strength of 



the second term is much greater than that of the first one (Quillen, 2003). 
Therefore, it is natural to perform a time-dependent shift <fi = ip — vt + 7 of 
the origin of the coordinate system. This shift makes the second resonance 
explicitly the guiding one. The resulting Hamiltonian is 

H = Af + (AS + v)j + Aej 1/2 cos V + A(3j 1/2 cos(ip - vt + 7). (16) 
Then we introduce the parameter A = AS + v and make a constant shift 



j = p — A/ (2 A), ip — cp — 7, reducing equation (16) to the form 

/ A\ 1/2 / A\ 1/2 

H = Ap i + Ae[p-—j caa<p + Ap[p- — ) cos(<p-ut), (17) 

and expand the coefficients in the neighborhood of p = 0, leaving only the 
lowest order (constant) terms. This gives 



H = Ap +Ae[ — 



2AJ 



cos <p + Ap I — 



2A 



cos(<£> — vt). 



(18) 
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Thus we have reduced the perturbed "second fundamental model for res- 
onance" (called so by Henrard & Lemaitre 1983), given by the Hamilto- 
nian ( JlSj ), to the perturbed "first fundamental model for resonance", given 
by the Hamiltonian (18), i.e., to the perturbed pendulum model. This has 
turned out to be possible because the guiding resonance (that emerging due 
to the spiral arms) is bifurcated in all cases, the librational "crescent" , born 
from the bifurcation, being situated always quite far from the origin of the co- 
ordinate system; see the phase space sections in figures 2-4 in Quillen (2003), 
figure 3 in |Quillen| ( |2009l ), and our Figures [I] and [2] The phase space sections 
in Figures [T] and [2] have been constructed by numerical integration of the 
eq uations of moti on with the Hamiltonian (fl~5|) in the same way as described 



in 



Quillen (2003); the coordinates in the sections are x 



' cos <t> and 



sin i 



In Figure [TJ the phase space section for model 6 of group C in | Quillen 



(2003) is shown. It corresponds to figure 4 (f) in Quillen (2003), except 
that we take here a much denser grid of initial data for the trajectories. In 
the case of Figure [TJ the value of 7 is equal to tc/2, as in Quillen (2003), 
while in the case of Figure [2] the value of 7 is chosen to be equal to zero. 
In both figures, the central regular island (small crescent) corresponds to 
the bar's resonance, while the outer regular island (big crescent) corresponds 
to the spiral's resonance. Setting 7 = results mostly in a shift in the 
angular location of the spiral's resonance in 0, that is why the relative angular 
orientation of the regular islands changes. The angular orientation of the 
center of the outer island (the fixed point of the spiral's resonance) changes 
by 7r/2. Apart from the change of the relative angular orientation of the 
regular islands, the phase space structure is qualitatively one and the same 
in Figures. The value of 7 does not influence the Lyapunov and diffusion 
timescales, because this parameter can be removed by a simple canonical 
transformation, namely, a constant shift in time. 

The chaotic domain is pronounced in both Figures. It emerges due to 
interaction of the bar's and spiral's resonances. The estimates of the Lya- 
punov and diffusion times are made in what follows just for such domains, 
emerging due to interaction of the resonances, in various model sets. 

Equation (18) describes the resonance duad, and, for estimating the Lya- 
punov time of the motion in the chaotic layer, we can apply formulas (10) 
and ( 14 ) in the cases of fast and slow chaos, respectively. 
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4 Lyapunov time estimates 



Comparing the Hamiltonians ( 1 ) and 



v, a 



\A\P(- 



V 

2A 



1/5 



18), one has: J 7 = —Ae 



b = ;Uo = lAWe] 1 / 2 (-25 



2A, 
A 



1/2 

i 

1/4 



g = 2A, n 

A = j^j, and e = — ~. The results of calculation of the Lyapunov times by 
formula (10) (when A > 1/2, i.e., in the majority of cases) and by formula 
(14) (when A < 1/2, only in two cases present) are given in Tables 1, 2, and 
3 for model groups A, B, and C, respectively. The model groups A, B, and 
C correspond to the sets of parameter values considered for the construction 
of the phase space sections in figures 2-4 in Quillen (2003). 

The Hamiltonian (15) has five parameters: 5, (3, e, v, and 7. All of 
them are expressed by means of formulas given in Quillen (2003) through 
original dynamical characteristics of the Galaxy and the solar neighborhood, 
such as the pattern speeds of the bar and spiral structure, the radius of the 
solar circle from the Galactic center, the radius at which the bar ends, etc. 
These five parameters elude straightforward interpretation in terms of the 
original dynamical characteristics, but they have straightforward meaning 
in the framework of the dynamical model (15): the parameters (3 and e 
characterize the strength of the bar's and spiral's resonances, respectively; 
5 is the "bifurcation control" parameter, as discussed in Quillen (2003); v 
measures the separation of the bar's and spiral's resonances in the frequency 
space; 7 is a phase shift, rather unimportant from the dynamical viewpoint, 
as discussed at the end of Section HI 

The models 1-10 of each model group correspond to the panel of 10 
Poincare sections in a corresponding figure in Quillen (2003): group A cor- 
responds to figure 2, group B to figure 3, and group C to figure 4. The 
parameters (3 and e are constants inside each model group. Their values are 
given in the notes to the Tables. The parameters 5 and v vary inside each 
model group. 

The obtained values of the Lyapunov time are all (except in one case, 
corresponding to adiabatic chaos) in the range from ^ 40 to ~ 80 time 
units. Inspection of the data given in Tables 1-3 also makes evident that the 
heuristic estimate ~ 2-njv (Holman & Murray, 1996 Quillen, 2003) severely 
underestimates the real Lyapunov time (by an order of magnitude in many 
cases). Besides, it is rather seldom that the considered dynamical systems, 
modeling the dynamics in the solar neighborhood, can be called "fully over- 
lapped" : the values of A (playing the role of the resonance overlap parameter, 
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see Section [2]) are typically greater than 1/2; so, chaos is non-adiabatic. 

Only model 9 (and, marginally, model 8) of model group B can be called 
adiabatic. The values of the problem parameters e, /3, u, and 5 in this 
case are such that the value of the separatrix map parameter A is less than 
1/2. Generally, the adiabaticity of chaos in any model can be checked by 
calculating A = \u\/uq, given the values of e, (3, u, and 5. 

As soon as our time unit is equal to one Galactic year divided by 2tt, 
we see that the obtained typical Lyapunov times are approximately equal 
to 10 Galactic years. How much is it, in comparison with, e.g., Lyapunov 
times of the solar system bodies, in comparable time units? The usual Lya- 
punov times for asteroids in the main belt are 3000-10000 years or more 
(Shevchenko, 2007), i.e., they are ~ 1000 asteroidal orbital periods or more, 
while the usual Lyapunov times for highly eccentric comets are of the order 
of one cometary orbital period (Shevchenko, 2007). So, the Lyapunov times 
of the chaotic motion in the considered Galactic model are much less (by at 
least two orders of magnitude) than the Lyapunov times for the main belt 
asteroids, but an order of magnitude greater than for the comets, if expressed 
in adequate time units. In general terms of the loss of predictability of the 
motion, the Galactic dynamical chaos is rather strong. 

From inspection of Tables 1-3 one can deduce that the Lyapunov time 
depends on the model parameters rather weakly being almost of the same 
order in all the models. Thus one can expect that choosing different values 
for the model parameters, such as relative strengths of the bar and spiral 
structures (e.g., choosing the spiral structure to be weaker), would not qual- 
itatively change the typical Lyapunov time value. However, note that the 
change of the model may radically reduce the extent of the chaotic domain. 
Then chaos is not "global" and so the probability that the Sun belongs to 
the chaotic domain of the phase space might be small. 

It is interesting that, as follows from the above estimates of Tl, the age of 
the Milky Way measured in its Lyapunov times is about 5-10. This means 
that now it is already practically impossible to restore exact initial conditions 
for the stellar dynamics in the solar neighborhood from any observational 
data. 
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5 Diffusion rates and ballistic flights 



Let us estimate the characteristic times of chaotic transport (called the diffu- 
sion times, when the diffusive approximation is used) in the chaotic domain 



of the phase space in the same model set of Quillen (2003). Knowledge of 



these times will allow one to judge on the possibility for migration of the 
Sun from the inner regions of the Milky Way to its current location. Such a 
possibility, arising due to overlapping of resonances in the phase space, was 



advocated and studied in detail by Minchev & Famaey (2010) and Minchev 



et al. (2011). To estimate the typical diffusion times, we shall base on the 



approach developed initially by Chirikov & Vecheslavov (1986, 1989) for the 



purposes of studies in cometary dynamics. Chirikov (1999) employed a sim- 
ilar approach in a study of the separatrix map dynamics (see page 11 and, 



in particular, formulas (20) and (21) in Chirikov 1999). 

First of all, a reservation should be made that it is only very approximate 
that we can consider the chaotic transport in the problem under study as 
diffusive. The matter is that the value of A, as follows from Tables 1-3 
the majority of models is rather low 
19991 p. 12-13) 



m 



A ~ 1. As explained in (Chirikov 



when A ~ 1, ". . .the layer width is reduced down to the 
size of a single kick. . . . Hence, the diffusion approximation becomes inap- 
plicable. Instead, the so-called ballistic relaxation comes into play which is 
much quicker. In other words, a slow diffusive motion ... is replaced now by 
rapid jumps of a trajectory over the whole layer . . .". (A "single kick" is 
the energy increment per iteration of the separatrix map.) Therefore all the 
diffusion rate estimates that we make in this Section should be regarded as 
extrapolation of diffusive description. In reality, they represent upper bounds 
for the real characteristic times of chaotic transport. 

According to Chirikov & Vecheslavov (1986, 1989), the diffusion rate (in 



the energy variable) in the main chaotic layer in the phase space of the Kepler 
mapQ approximately equals the mean (over time) squared energy increment 
per iteration of this map. Analogously, in the case of the ordinary separatrix 
map ([2]), the diffusion rate (in the energy variable) in the main chaotic layer 
in the phase space of the map approximately equals the mean squared energy 
increment, i.e., {W 2 sin 2 Tj). Averaging over the interval < r, < 2n, we find 



the diffusion rate to be D 



map 



W 2 /2. 



1 The Kepler map is a kind of a general separatrix map, the parabolic motion playing 



the role of a separatrix (Shevchenko 2010). 
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In the case of the separatrix algorithmic map, the chaotic layer com- 
ponents corresponding to prograde rotations, retrograde rotations, and li- 
b-rations of the phase variable should be considered separately. In the two 
(prograde and retrograde) rotation cases the diffusion rate -D map in the en- 
ergy in the map phase space obviously equals ~ (W + ) 2 /2 and « (W~) 2 /2 for 
the prograde and retrograde rotations, respectively. Employing the formulas 
for n, uq, a, b (b — 0), A, and e, given at the beginning of Section [ij one can 
calculate the parameters W + and W~ of the separatrix algorithmic map (j5|. 
If A > 1/2, the equality 6 = implies \W~\ <C \W + \, while a = implies 
\W~\ S> The component of the chaotic layer corresponding to reverse 

(or direct) rotations does not exist, if W~ (or, respectively, W + ) is equal 
to zero; its measure is zero. The other circulation component with non-zero 
measure is described by the ordinary separatrix map ^ with the parameters 
A and W = W (the non-zero value among W + and W~); its extent in w 
is ~ A|W| (the half-width of the chaotic layer in the case of fast chaos; see 



Shevchenko 2008a). 



Consider the libration component of the chaotic layer. Then W~ and W + 
alternate (replace each other) at each iteration of the separatrix algorithmic 
map ([5]). It is straightforward to show (Shevchenko, 2007) that, if W~ or W + 



is equal to zero, the separatrix algorithmic map (|5j) on the doubled iteration 
step reduces to the ordinary separatrix map ^ with the doubled value of A 
and the value of W equal to a non-zero value of W ± . Taking into account 
that one iteration of the new map corresponds to two iterations of the original 
one, the diffusion rate referred to the original map time units is 

Anap^(W^) 2 , (19) 

where W is the non-zero value among W + and W~ . 

For the circulation component of the layer, one has D map ~ (W + ) 2 /2 (if 
b = 0) or Anap ^ (W-f/2 (if a = 0). 

As soon as the libration motion is reducible to the ordinary separatrix 
map (J2J) with the doubled value of A, the layer's extent in w on the side 
of librations doubles, it becomes ~ 2A|W|. Note that the parameters A 
and W are considered here as independent from each other. Therefore, the 
chaotic domain corresponding to libration dominates in extent, and for a 
rough estimate of the diffusion rate over the entire layer it is sufficient to 
make an estimate for the libration component alone. 
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In the case of the Hamiltonian (18) b = and 77 = 0, so, W = W~ 
where 



2 exp(vrA/2) 



W + = e\A 2 (X) = AttsX 2 V ( 20 ) 

smh(7rA) 



(see equation ([7])), therefore 

sinh(7rA) 



D^^ll^^^) (21) 



in the lib ration case. 

To obtain the diffusion rate referred to real time units, it is necessary 
to transform the map time units into the real ones. This is performed by 
dividing the diffusion rate referred to map time units by the mean period of 
phase rotations (half-periods of librations) inside the chaotic layer, because 
this mean period is nothing but the average time interval corresponding to 
one map iteration. Consequently, the diffusion rate referred to real time units 



is D = \n\D map /T sx , where T sx is given by formula (13). 

We define the characteristic diffusion time across the chaotic layer to be 
equal to the inverse of the diffusion rate. Therefore, it is just the time needed 
for the diffusing particle to cover the relative energy interval equal to one. 
Note that the maximum possible deviation in the relative energy w from zero 



in the libration case is equal to —2 (Chirikov, 1979); therefore, the defined 



diffusion time gives an appropriate time estimate for the global mixing inside 
the chaotic layer, of course, if the chaotic layer is broad enough. 



In our Hamiltonian (18) b = 0, so, W ± = W + , and one gets for the 
diffusion time 

1_ T SX (A,W/+) ^ 4T SX (A,H/+) 

d D ~ |Q|Anap " M(W + ) 2 ' 1 ' 



where 



32<= 

T sx (A,W/+)^Aln— — (23) 



see equation (13)), e is the base of natural logarithms, and W + is given by 



formula (20). Finally one has 



4A , 32e , „ 

T ^mw^ ln w^\- (24) 
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In the case of the prograde rotation component of the chaotic layer, the 
diffusion rate D map ~ (W + ) 2 /2; therefore, the diffusion time is two times less 
than that given by formula (24). 

The results of the calculation of the diffusion times by formula (24) 
are given in Tables 4, 5, and 6 for model groups A, B, and C, respectively 

Inspection of the data in Tables 4-6 makes it evident that the diffusion 
times vary considerably inside the model groups: by as much as four orders of 
magnitude in group C. Taking into account that our time unit is equal to one 
Galactic year divided by 2tt, i.e., our time unit ~ 32 Myr (1 Myr = 10 6 year; 
1 Galactic year 200 Myr), it is clear that the obtained diffusion times in 
most of the models with ordinal numbers up to three are greater than 10 Gyr 
(1 gigayear = 1 Gyr = 10 9 year), making large-scale radial chaotic migration 
improbable. Such "junior" models all have positive values of 5. However, 
models 7-9 in group B also have large values of 



Minchev & Famaey (2010) found in detailed numerical experiments that, 
due to the resonance overlap of the bar and spiral structure, the Galactic 
disk mixes in about 3 Gyr. From our Tables 4-6 it is obvious (taking into 
account that our time unit ~ 32 Myr) that most of the models with negative 
values of 5 provide chaotic mixing effective enough to be consistent with 
the numerical-experimental results of Minchev & Famaey (2010). Almost all 
models with large ordinal numbers, except models 7-9 in group B, permit 
such migration. Generally, as follows from data in Tables 4-6, the diffusion 
rate strongly depends on the radial position in the Galaxy. 

Note that the Hamiltonian model (15) assumes that the guiding radius 
in the unperturbed dynamics is fixed. If migration is present, the guiding 
radius varies; so, a more refined model should be used to account for this 
variation self-consistently. One should also mention that chaotic transport 
can be ubiquitous in the Galaxy, though its dynamical origin is presumably 
different from that considered here. Quillen et al. (2010) showed that the 



diffusion might occur in many regions due to interaction of multiple waves, 
and so overlapping of resonances should be common; besides, if the bar slows 
down (e.g., Weinberg & Katz 2007 and references therein) resonances can 
sweep through the Galaxy (Quillen et al. 2010). Of course, these processes 
cannot be described in the specific dynamical model considered here, but this 
model already provides an insight into the possible effectiveness of chaotic 
transport in the Galaxy. 

The most important conclusion of this Section is that models that permit 
large-scale radial chaotic migration of the Sun (from the inner regions of the 
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Milky Way to its current location) do exist. This confirms the dynamical 
possibility of the migration concept advocated by Minchev & Famaey (2010) 
and Minchev et al. (2011). What is more, due to the possibility of ballistic 
flights mentioned in the beginning of this Section, the chaotic mixing might 
be far more effective and quicker than in the case of normal diffusion. Ob- 
viously, the effect of ballistic relaxation should be explored in detail in the 
future. 



6 Discussion and conclusions 



We have considered how the Lyapunov and diffusion timescales of the stel- 
lar dynamics in the solar neighborhood can be estimated. We have used 
Quillen's (2003) model to describe interaction of the "spiral" and "bar" non- 
linear resonances in the phase space of the motion. A method of analytical 
estimation of the maximum Lyapunov exponents of the orbital motion has 
been applied to the solar neighborhood dynamics. The analytical treat- 
ment has been performed within a framework of the separatrix map theory 



(Shevchenko, 2000a, 2002, 2007), describing the motion near the separatrices 



of a perturbed nonlinear resonance. 

The Lyapunov times turn out to be basically in the range from 6 to 13 
Galactic years. In comparison with the Lyapunov times of the solar system 
bodies (made in adequate time units), the Galactic dynamical chaos is rather 
strong in general terms of the loss of predictability of the motion. An in- 
teresting inference is that, as soon as the age of the Milky Way measured 
in its Lyapunov times is about 5-10, now it is already practically impossi- 
ble to restore exact initial conditions for the stellar dynamics in the solar 
neighborhood from any observational data. 

We have estimated also the diffusion times, based on the approach de- 



veloped initially by Chirikov & Vecheslavov ( 1986 , 1989 ) for the purposes of 



studies in cometary dynamics. We have found that, in a number of mod- 
els, the diffusion times turn out to be small enough to permit radial chaotic 
migration of the Sun from the inner regions of the Milky Way to its cur- 
rent location. In other words, dynamically adequate models that permit 
large-scale radial chaotic migration do exist. This confirms the dynamical 



possibility of the migration concept advocated by Minchev & Famaey (2010) 



and Minchev et al. (2011). Due to the possibility of ballistic flights inside 
the chaotic layer, arising because A ~ 1, the chaotic mixing might be even 
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far more effective and quicker than in the case of normal diffusion. 

We have shown that only in a narrow range of possible values of the 
problem parameters e, /3, u, and S the Galactic chaos is adiabatic because the 
values of the separatrix map parameter A, playing the role of the resonance 
overlap parameter, are typically greater than 1/2; in other words, adiabatic 
chaos (A < 1/2) seems to be not characteristic for the dynamics in the solar 
neighborhood. 

The author is thankful to Alice Quillen for advice and comments. It is 
also a pleasure to thank the referee for helpful remarks. 
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Table 1: Lyapunov time estimates for model group A (e = —0.004, (3 = 
0.0006; e = 0.15; this model group corresponds to the graph panel in figure 2 
in |Quillen| (|2003|)) 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


S 


0.068 


0.034 


0.018 


0.001 


-0.006 


-0.009 


-0.016 


-0.019 


-0.032 


-0.049 


V 


1.000 


0.750 


0.625 


0.500 


0.450 


0.425 


0.375 


0.350 


0.250 


0.125 


A 


4.07 


3.13 


2.65 


2.15 


1.94 


1.84 


1.64 


1.53 


1.11 


0.56 


T L 


44.3 


39.7 


38.5 


38.3 


38.8 


39.3 


40.5 


41.5 


48.2 


74.2 



Table 2: Lyapunov time estimates for model group B (e = —0.004, (5 = 
0.0005; e = 0.125; this model group corresponds to the graph panel in figure 3 
in 



Quillen (2003)) 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


S 


0.068 


0.034 


0.018 


0.001 


-0.006 


-0.009 


-0.016 


-0.019 


-0.039 


-0.066 


V 


0.700 


0.465 


0.347 


0.230 


0.183 


0.159 


0.112 


0.089 


-0.052 


-0.240 


A 


3.37 


2.32 


1.78 


1.20 


0.97 


0.85 


0.60 


0.48 


0.29 


1.42 


T L 


49.4 


46.8 


49.0 


57.1 


64.3 


69.9 


87.5 


69.6 


123.7 


60.4 



Table 3: 
0.0006; £ 
in 



Lyapunov time estimates for model group C (e = —0.002, (3 = 
= 0.3; this model group corresponds to the graph panel in figure 4 



Quillen (2003)) 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


S 


0.068 


0.034 


0.018 


0.001 


-0.006 


-0.009 


-0.016 


-0.019 


-0.032 


-0.066 


V 


1.200 


0.840 


0.660 


0.480 


0.408 


0.372 


0.300 


0.264 


0.120 


-0.240 


A 


6.44 


4.78 


3.89 


2.95 


2.55 


2.35 


1.93 


1.72 


0.82 


2.01 


Tl 


76.4 


60.4 


53.2 


47.4 


46.1 


46.0 


47.0 


48.6 


74.1 


60.9 
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Table 4: Diffusion time estimates for model group A 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


A 


4.07 


3.13 


2.65 


2.15 


1.94 


1.84 


1.64 


1.53 


1.11 


0.56 


W + 


0.104 


0.271 


0.412 


0.595 


0.673 


0.709 


0.771 


0.798 


0.813 


0.506 


T d 


7900 


1100 


440 


210 


160 


140 


120 


120 


120 


400 



Table 5: Diffusion time estimates for model group B 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


A 


3.37 


2.32 


1.78 


1.20 


0.97 


0.85 


0.60 


0.48 


0.29 


1.42 


W+ 


0.179 


0.442 


0.608 


0.687 


0.646 


0.600 


0.451 


0.358 


0.200 


0.681 


T d 


3000 


450 


240 


210 


250 


310 


610 


1000 


4100 


230 



Table 6: Diffusion time estimates for model group C 



Model 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


A 


6.44 


4.78 


3.89 


2.95 


2.55 


2.35 


1.93 


1.72 


0.82 


2.01 


W+ 


0.013 


0.094 


0.253 


0.638 


0.893 


1.04 


1.35 


1.50 


1.41 


1.30 


T d 


9.4 x 10 5 


13000 


1600 


230 


110 


84 


49 


41 


60 


70 
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Figure 1: Typical example of the phase space section (model 6 of group C). 
The inner regular island (small crescent) corresponds to the bar's resonance 
and the outer regular island (big crescent) corresponds to the spiral's reso- 
nance. 
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Figure 2: Same as in Figure [TJ but 7 = 0. Apart from the change of the 
relative angular orientation of the regular islands, the phase space structure 
qualitatively remains the same. 
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